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Abstract 

We report variational and diffusion quantum Monte Carlo (VMC and DMC) studies of the 
binding curve of the ground-state chromium dimer. We employed various single determinant 
(SD) or multi-determinant (MD) wavefunctions multiplied by a Jastrow fuctor as a trial/guiding 
wavefunction. The molecular orbitals (MOs) in the SD were calculated using restricted or unre- 
stricted Hartree-Fock or density functional theory (DFT) calculations where five commonly-used 
local (SVWN5), semi-local (PW91PW91 and BLYP), and hybrid (B1LYP and B3LYP) functionals 
were examined. The MD expansions were obtained from the complete-active space SCF, gener- 
alized valence bond, and unrestricted configuration interaction methods. We also adopted the 
UB3LYP-MOs to construct the MD expansion (UB3LYP-MD) and optimized their coefficients at 
the VMC level. In addition to the wavefunction dependence, we investigated the time-step bias 
in the DMC calculation and the effects of pseudopotentials and backflow transformation for the 
UB3LYP-SD case. Some of the VMC binding curves show a flat or quite shallow well bottom, 
which gets recovered deeper by DMC. All the DMC binding curves have a minimum indicating a 
bound state, but the comparison of atomic and molecular energies gives rise to a negative binding 
energy for all the DMC as well as VMC calculations. 
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I. INTRODUCTION 



The chromium dimer (Cr 2 ) has attracted a lot of attention as a prototype to understand 
the d-d binding in both experimental jl-7] and theoretical studies [9-50]. The ground state 
is experimentally observed to be a singlet state, whereas the ground state of the 

constituent Cr atom 3d 5 As 1 ( 7 S) has the highest spin multiplicity in the 3c? atoms. Recent 
spectroscopic experiments reported an equilibrium bond length (R e ) and binding energy 
(D e ) of 1.6788 A|2j and 1.53 ± 0.05 eV[3 , respectively (some older experiments reported R e 
of 1.68 AQ and D e of 1.44 ± 0.05 eV j4] and 1.42 ± 0.10 eV However, no theoretical 
study has provided a quantitatively satisfactory result for Cr 2 so far, since its chemical 
binding is rather complicated. 

From the theoretical viewpoint, there are two extremely different pictures to understand 
the chemical binding in Cr 2 qualitatively. The first one is based on elementary molecular 
orbital (MO) theories. In this framework, Cr 2 is treated as a closed shell (single-determinant) 
state, with all the bonding orbitals occupied (\o 2 g '2o 2 g lir^lSg, arising from the 3c? and 4s 
orbitals). Cr 2 is therefore interpreted as having a sextuple bond, which is the highest 
multiple bond in any diatomic molecule. This is naively consistent with the observed short 
bond length (w 1.68 A) , whereas the experimental D e (« 1.5 eV) is rather small in 
the context of multiple bonds, which is smaller than that of singly bonded Cu 2 (~ 2.0 eV). 
This fact may imply that an elementary one-electron approximation picture is invalid for 
Cr 2 and hence "non-dynamical" (static) correlation effects are important. Indeed, restricted 
Hartree-Fock (RHF) does not only give an incorrect dissociation behavior, but also gives 
rise to a ground state energy of Cr 2 far above that of the two isolated atoms (about 20 eV 
higher) 13]. 

The second picture, which is the extreme opposite to the first one, emphasizes the relative 



stability of high-spin atomic states 



13] . This is the molecular analogue of antiferromagnetism 



and can be treated at the crudest level of theory by the unrestricted (or broken symmetry) 
Hartree-Fock (UHF) theory using a spin unrestricted single-determinant of symmetry broken 



molecular spin orbitals [30|, |48|. The method can approximately deal with the static correla- 
tion effects. Up and down-spin electrons are antiferromagnetically localized on each of the Cr 
atoms, where the charge density has a symmetry but the spin density has a Cqo V sym- 
metry. The wavefunction represents the nearest possible single-determinant approximation 
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to the X E+ state, but it is not an eigenfunction of total spin, leading to spin contamina- 
tion. This can be easily remedied by a properly symmetrized multi-determinant expansion 
of the UHF wavefunction, i.e., the generalized valence bond (GVB) method. However, the 
GVB method gives rise to a very weakly bound molecule (D e = 0.35 eV) with a very long 
bond length (R e = 3.1 A) 13f]. Although a complete active space (CAS) self-consistent field 
(SCF) method can also consider the static correlation, a typical CASSCF with a 12 orbital 
active space and 12 active electrons, CASSCF(12,12), provides a poor result of R e and D e , 
similar to GVB [48]. These disagreements with experiments indicate dynamical correlation 
effects are also important for achieving quantitatively satisfactory results. In summary, the 
chemical binding in Cr 2 involves a highly complicated blend of 4s-4s and 3d-3d interactions 
with antiferromagnetic coupling. 

To understand such a complicated binding mechanism, a large number of ab initio studies 
have been performed based on either traditional MO theories 9H32I] or density functional 
theories (DFT) 33m0j]. In MO, both single- and multi-reference theories were studied. 
Within the single-reference theory, the coupled-cluster approach with single and double 
excitations including triples noniteratively, CCSD(T), is one of the most accurate methods. 
Restricted CCSD(T) calculations give a very weak binding (D e = 0.38 eV) with a short bond 
length (R e = 1.60 A), whereas unrestricted CCSD(T) calculations provide a better binding 
energy (D e = 0.89 eV), but with a longer bond length (R e = 2.54 A). In multi-reference 
theories, a typical reference is a CASSCF(12,12) calculation which is responsible for the 
static correlation. Using the CAS reference, dynamical correlation can be taken into account 
by multi-reference configuration interaction (MRCI) 18| . coupled cluster (MRCC) 2Jj, or 
second-order perturbation (CASPT2) methods (23], Although these multi-reference methods 
give a better R e than single-reference methods, there is still room to improve further the 
accuracy of D e . 

In DFT, various exchange-correlation (XC) functionals are available such as the local 
(spin) density approximation (LDA/LSDA), the generalized gradient approximation (GGA), 
as well as the hybrid XC functionals. Both restricted LDA (RLDA) and the unrestricted 



formalism (ULDA) give rise to overbinding, i.e., too short R e and too large D e 



33-37 



is a well-known failure of LDA. Both restricted and unrestricted GGA calculations 



^jyhich 



38 



47l | generally improve the LDA discrepancy, but it is difficult to choose the "best" GGA 



functional because one may give a better R e , but it may give a worse D e , and vice versa. 
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Restricted B3LYP, though it is popular for covalent molecules, gives rise to an unbound 



molecular state 



38 



48 



Although unrestricted B3LYP reproduces a bound state, it provides 

. These results may imply 
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a smaller D e of ~ 4.0 eV with much larger R e of ~ 2.5 A 
that the difficulties for the binding of Cr 2 would originate from a delicate balance between 
exchange and correlation in DFT for chemical binding of Cr 2 . 
Quantum Monte Carlo (QMC) methods 
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52j are one of the most accurate techniques 



in state-of-the-art ab-initio calculations for quantitative descriptions of electronic structures. 
There are two typical QMC calculations, i.e., variational and diffusion Monte Carlo (VMC 
and DMC) methods. VMC is not usually accurate enough since its result strongly depends 
on the correlated trial wavefunction adopted. DMC is a technique for numerically solving 
the many- electron Schodinger equation for stationary states using imaginary time evolution. 
The fixed-node approximation is usually assumed to maintain the fermionic anti-symmetry 
in DMC. Although the fixed-node DMC can accurately evaluate the ground state energy 
of many atoms and molecules using only the trial node from a single determinant (SD), it 
sometimes fails, especially for "near-degeneracy" systems such as the Be atom. This implies 
that the fixed-node DMC method can work well for the dynamic correlation, but not for 
the static correlation which should be included at the stage of choosing the fixed-node trial 
wavefunction. The Cr 2 molecule is also considered as a near-degeneracy molecular system 
and hence a good challenge to QMC. 

In this study we performed VMC and DMC calculations of Cr 2 with several choice of trial 
wavef unctions. Variety of the XC functionals are examined to construct orbital functions 
including HF, SVWN LDA, PW91 and BLYP GGA, B1LYP and B3LYP hybrid functionals, 
with both restricted and unrestricted treatments. In addition to a single-determinant form 
of the many-body wavefunction, we also tried several multi-determinant (MD) forms. For 



the orbital part we also introduced the backflow transformation 



66 



68 



. Performance of the 



XC functionals are examined based on the variational principle with respect to the fixed 



node 



51 



52|. The chemical binding of the ground-state Cr 2 may be examined in two ways: 



(i) use of the binding curve and (ii) a comparison of atomic and molecular energies. 

The present paper is organized as follows: Section II describes the present computational 
methods. Section III provides numerical results and discussion. Section IV summarizes the 
present study. 
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II. COMPUTATIONAL METHODS 



We replaced the inner Neon core (10 core electrons) of the Cr atom with a small core 
norm-conserving pseudopotential which is constructed from Dirac-Fock atomic solutions. 
More specifically we employed the Lee-Needs (LN) soft pseudopotential 56J using Troullier- 
Martins construction. To check the dependence on the potential we also used a small core 
Burkatzki-Filippi-Dolg (BFD) pseudopotential {stJ. An electronic many-body wavefunction 
is composed of anti-symmetrized products of orbital functions which is generated from DFT 



or HF. Various XC 
the nodal surface 



unctionals give different orbitals from which different structures of 
53| of the many-body wavefunction (the surface in configuration space 
on which the wavefunction is zero and across which it changes sign) are constructed. For 
different nodal structures, we can exploit the variational principle to see which choice is 
better, comparing the ground state energy 54j. In this study we tested several combinations 
of (i) the choice of XC, (ii) spin restricted/unrestricted treatment of orbitals, (iii) SD/MD, 
(iv) choice of pseudopotentials, and (v) with/without backflow degrees of freedom 66|. 

In VMC the ground state energy is evaluated as the expectation value of the Hamiltonian 
H with a many-body trial wavefunction, 

jm*Hmd& _ J \y\ 2 y- 1 HydR 

where R = (ri, tn) denotes an electronic configuration of valence electrons in a molecule 
(N = 28 with 14 up and down spins, respectively, for Cr2), and the energy has been written 
as an average of the local energy, El (R) = ^~ l H^> , over the probability distribution p(R) = 
\^\ 2 / J |\I/| 2 <iR. The energy expectation value is evaluated from Monte Carlo integration, 
using the Metropolis algorithm to generate electronic configurations distributed according 
top(R). 

In the DMC method the ground-state component of a trial wavefunction which overlaps 
with the exact one is projected out by evolving an ensemble of electronic configurations using 
the imaginary-time Schrodinger equation. Attempts to carry out this procedure exactly 
result in a "fermion sign problem" , which is circumvented by constraining the nodal surface 
of the wavefunction to equal that of the trial wavefunction. The DMC energy calculated 
with this fixed-node constraint is equal or higher than the exact ground-state energy, and 
becomes equal to the exact one if and only if the fixed nodal surface is exact. 
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We used the many-body wavefunction taking the form, 



(R) = e J W ■ Fas (r 1 ,---,r N ) 



(2) 



where Fas (ri, • • • , rjv) is the anti-symmetrized products of orbital functions {ipf (ij)} such 
as a Slater determinant. In this study the function is expanded as 



F AS (ri, • • ■ , r^) = ^ cy ■ Dj (f l9 • ■ ■ , r#/ 2 ) £>j (?iv/2+i, ■ ■ ■ , ?at) 



(3) 



i=0 j=0 



by determinants with coefficients c^-. Fq ''^ (fi, • • • , f jv/2) is the ground-state single Slater 



determinant formed only by the occupied orbitals for each spin, and -DJ_^ = F„ 



U) 



^0 



(i) 

corresponds to excited state configurations, where Pn-^m denotes the excitation from the 
occupied state ( r j) the virtual state ^(r,,). In Eq. (j3J) arguments with a tilde 



r i + 6 ( r i) " " " > r iv), denote the backflow shift |66| 



JV MD 



corresponds to SD, 



and d 



Ci ' Oij } D 



t 



Dj to a spin restricted wavefunction. The orbital functions, 



{ipl (ij-)}, were expanded with a contracted Gaussian basis set (17sl8pl5d6fV 



Gaussian03 



8s8p7d3f]. 

was used for SCF calculations, while we used CASINO ver.3.0 59] for QMC 
calculations. Some calculations (HF and GVB) were carried out using GAMESS 60[ for 
SCF and Q Walk |6l| for QMC. We attempted to use another many-body wavefunction form, 



Pfaffian 



741 ] . which is available in QWalk. However an optimization procedure for the off- 



diagonal elements of the Pfaffian did not work well and then we could not obtain reliable 
results, so not reported here in detai . 



As for the Jastrow factor e J( - n > 55|, the function J (R) is given as 62 ] 



A' 



N-l N N ions 
i=l j=i+l 1=1 i=l 



N ionB N-l N 
1=1 i=l j=i+l 



(4) 



where suffices i,j and I specify electronic and ionic positions, respectively. Each term, u, 
X, and /, takes into account the dynamical correlation due to electron-electron, electron- 
nucleus, and electron-electron-nucleus coalescence, respectively. Cutoff lengths are intro- 
duced to make each term quadratically fall off to zero at the radius. Specific forms used 
in each QMC code are described in the appendix. The electron-electron cusp condition 
63j is imposed in the it-term. Variational parameters in J (R) are designed to be able to 
include spin polarized case and are optimized by VMC individually at each bond length 
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with fixed cutoff lengths. Those parameters were optimized by the variance minimiza- 



tion 64( as well as the energy minimization 75] procedures. The backflow shift 66], 
6 ({ r i}) = t,f e ({ r ij}) +^ Ar ({ r i/}); is introduced to modify the nodal surface variation- 
ally. Each term, ee and eN, is expanded with the power of inter-particle distances (up to 
8th order in this study) with proper cutoff radii which are fixed in the same values as those 



67 



in u and \ in our Jastrow factor 
by the filtered reweighted variance minimization scheme 
of freedom. 



681 ] . The parameters in the backflow were optimized 



64] allowing spin polarized degrees 



For the singlet Cr 2 ground state[l-7] the spin restricted SCF treatment seems not to 



give a proper description of localized spin polarization on atomic sites [69] . Even using the 
restricted nodal surface, however, the DMC projection can still recover the proper localized 
spin polarization. Indeed we found the RHF nodal surface gives a variationally better 
description than the UHF one. We therefore investigated spin restricted methods as a 
possibility as well as unrestricted ones. In order to generate the nodal surface, we employed 
five commonly-used XC functionals, SVWN5(LDA), PW91PW91 and BLYP (GGA), B1LYP 
and B3LYP (hybrid), shown in Table. [H 

For DMC we started statistical accumulations of 2,000 walkers after equilibration of 1,000 
steps with 5t= 0.01 [a. u.], for which we have confirmed that the time step bias is within the 
statistical noise considered here. The present non-local pseudopotentials were evaluated by 
the T-move scheme [70J which is devised to reduce the instability and bias due to the locality 
approximation 7^]. 



III. RESULTS AND DISCUSSION 



A. Single determinant calculations 



Binding curves obtained from SCF and QMC calculations using unrestricted SD wave- 
functions are shown in Fig. □ (a) - (c). At the SCF level, LDA (USVWN5) and 
GGA(UPW91PW91 and UBLYP) recover a bond length, R e , near to the experimental 



value, 3.1748 [a.u.][46], while the other hybrid functionals give a longer R e . At the QMC 
level, using any of the XC functionals, however, the bond length, R e , is very similar and 
much overestimated compared to experiment. The results imply that the non-local HF ex- 
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change favors a longer R e . QMC essentially takes the non-local exchange into account even 
if the trial/guiding wavefunction is constructed without the HF exchange. VMC gives a 
quite shallow bottom of the binding curve, while DMC gets the curve recovered deeper. 



However, DMC sti 
D e , 0.0541 hartree 



1 unsatisfactorily recovers about 30% of the experimental binding energy 



46]. 



At the experimental R e some of the SCF calculations give lower molecular energies than 
twice the isolated atomic energies (zero-binding energies shown in Table IHIj) , as shown in 
Table M LDA(USVWN5) and GGA(UPW91PW91) recover 164.5 % and 70.3 % of the 
experimental D e , respectively, while UB3LYP gives 3.7 %. Underlined values in Table [TT] 
highlight those lower than the zero-binding energies, which means that the molecule is 
bound at the experimental R e . For example, RB3LYP-SCF does not bind the molecule, 



which is consistent with previous works 



38 



481 ] . To compare the atomic and molecular 



energies, the same basis sets and pseudopotentials were used. The RHF, UHF, and QMC 
calculations using the RHF and UHF trial wavefunctions were obtained by GAMESS/QWalk 
while the others are by Gaussian03/CASINO. A lower DFT-SCF energy does not mean a 
better description of the total energy because it does not necessarily satisfy the variational 
principle, while a lower QMC energy implies a solution closer to the exact one because of 
its variational property of the total energy. In Table [Til therefore, only the best (lowest) 
zero-binding energies of VMC and DMC using the B1LYP trial/guiding wavefunction are 
shown. Hence none of the QMC calculations at the SD level can reproduce a bound state 
at the experimental R e . 

It is worth noting in Table UTI that RHF gives a lower DMC energy than UHF. At the SCF 
level, in turn, RHF gives a much higher energy than any other, supporting the consensus 
that restricted treatments cannot well describe the localized amplitude of the wavefunction 
at ionic sites of a spin polarized system 69(. In DMC, however, the amplitude of the many- 
body wavefunction is automatically adjusted by the projection operation and is not directly 
governed by the XC approximation. The results show that the RHF nodal surface is superior 
to the UHF one, leading us to examine restricted methods for the generation of the trial 
nodal surface. 

Figure |2] (a) - (c) shows binding curves evaluated by restricted methods. At the SCF 
level, they well reproduce the experimental R e . A quite deep well bottom gives rise to an 
overestimate of D e , which is attributed to a well-known failure of restricted methods 69 ]. 



S 



Such an overestimate is modestly improved in DMC because the imaginary projection re- 
laxes the many-body wavefunction at larger distances. From the viewpoint of evaluating 
R e , RB3LYP gives the best restricted fixed nodes, though it is found to be variationally 
worse than any of unrestricted nodes. The variationally optimal fixed nodes within the SD 
approximation are obtained from UB3LYP, though it overestimates R e by around 50 % and 
underestimates D e by 30 %. The backflow transformation improves the ground state energy, 
but makes no significant improvements on R e and D e , as shown in Fig. [8] (a) and (b). 

We also investigated the dependence on pseudopotentials for the best fixed nodes within 
the SD approximation, i.e., we performed UB3LYP-QMC with backflow for two different 
potentials, LN and BFD. The latter uses a contracted Gaussian basis set (33s29pl9d2f)/ 
[5s5p4d2f] {57]. Though both pseudopotentials have the same effective valence electrons, 
the atomic energies are quite different from each other, as shown in Table [IV] in terms of 
zero-binding energies. We note that considerably faster SCF convergence is observed with 
BFD, by a factor of five faster than LN, while the QMC statistical qualities (statistical noise, 
auto-correlation, and population fluctuation) are almost the same for both potentials. The 
LN/UB3LYP-SCF calculation gives a weak binding at the experimental R e , but such a bound 
state disappears for BFD/UB3LYP-SCF (unbound), as seen in Table HVl Comparisons of 
binding curves are shown in Fig. [3] (a) - (c). There is no significant difference between 
LN and BFD in the predictions of R e and D e , though BFD slightly underestimates D e and 
overestimates R e . In summary, this result may justify our choice of the LN pseudopotential. 

The time step dependence of the DMC energies is shown in Fig. HI It is confirmed that 
the result with dt = 0.01 agrees with the results using other choices of dt to within one 
standard deviation a. The sudden decrease in the energy at dt = 0.001 is found to be 



similar to Fig. 7 of Ref. 



72|. 



B. Restricted multi determinant calculations 



Valence bond (VB) type many-body wavefunctions are expected to give a proper descrip- 



tion of the spin polarized Cr 2 



7lj with a compact form of the MD expansion. Generalized 



VB (GVB) SCF is available in GAMESS|60| and we can use GVB orbitals to generate QMC 
trial/guiding wavefunctions. Table PVl shows the symmetries of the UHF natural orbitals 
(NO) near the HOMO-LUMO level. We considered 12 active occupied molecular orbitals 
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up to level 20, which arises from the 4s(l)3d(5) atomic orbitals. The GVB function is formed 

as 



(6) _ a J $(8) . TT *(2) 
GVB — A \ ^corc 11 ^GVB;p 

p=l 



(5) 



(S) (2) 

where A denotes an antisymmetrizer and $core the core contribution. $gvb-p consists of the 
(GVB) orbital pairs, for which each of j-occupied orbitals (j G {9, • • • , 14}) in the active 
space is paired with a virtual orbital with the same symmetry. Hence the following six pairs 
are involved: p = (10, 19), (11, 18), (14, 15), (13, 16), (12, 17), and (9, 20). The orbitals in the 
active space were optimized by a SCF procedure with respect to ^^vb- An ex phcit form of 



^gvb takes a restricted CI expansion, D] = Dj, dj 



Ci ' ^iji aS, 



7 GVB 
AS 



i=0 



gg ■ 1 + eg • -Pg^Sg 
9io ■ 1 + eio " -Plo->5i 



<8> [gu ■ 1 + e u ■ Pi4^s 14 



(6) 



where gj and Cj are coefficients such that c = g 9 ■ g w ■ ■ ■ g 14 etc. Sj denotes the index of 
a virtual orbital with the same symmetry as the j-occupied orbital. The operators, 1 and 
Pj^Sj, are the identity and permutation operators, respectively. Pj^Sj swaps the j-occupied 
orbital in Dq into the 5j-virtual one. As discussed later a usual CI treatment in a quantum 
chemistry code gives hundreds of thousands of terms in the expansion, all of which can not 
be included in a QMC calculation. GVB provides, in contrast, a very compact form of the 
MD expansion with only 64 terms. 

Starting with UHF-NO and optimizing the coefficients, gj and e^, as well as the orbitals, 
GVB-SCF gives a variationally better description than RHF within a restricted SCF treat- 
ment, as seen in Table ED With the coefficients given by GVB-SCF, the wavefunction 
achieves a better (lower) energy than HF at the VMC level, but it turns out to be higher 
than HF using DMC (the row of 'GVB (6)' in Table I VII) . Then we tried to optimize the 
coefficients further by VMC. A total of 64 coefficients can be reduced to 48 independent 
variables byits symmetry. We adopted a mixed scheme between energy and variance mini- 



mization 



75[ with 95% weight on the former. Even though we ignore the above symmetry 



reduction to optimize 64 parameters independently, the optimized values of the parameters 
roughly satisfy their symmetries. Using the GVB nodes with these coefficients, we obtained 
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a better DMC value than when using the HF nodes, but still above the zero-binding energy 
(the row of 'GVB(6)opt' in Table EI}. 

Next we considered a restricted CASSCF node (I)t = Dj, = Cj, ■ Sij), having the form 

of 

^CAS-l AfeAS-l 

^a C s AS = E «"Df= E ^^o- (7) 

i=0 i=0 

An initial guess for the orbitals was taken from a UHF-NO calculation and optimized self- 
consistently with respect to the above many-body wavefunction. We employed (restricted) 
CASSCF(2,4) and CASSCF(2,7) methods, in which the number of expansion terms amounts 
to A r cAs=16 and 49, respectively. The results are shown in Fig. [5] (a) - (c). Though 
they could not achieve a variationally lower energy than UB3LYP-SD in terms of the final 
DMC energy, several interesting behaviors are found as follows: CASSCF(2,7)-SCF gives 
a binding curve with a similar shape to UB3LYP-SCF, though overestimating R e . At the 
QMC level, in turn, the evaluated value of R e gets shorter, and the shape of the binding 
curve is similar to the restricted SD cases. This implies that the terms in the MD expansion 
well describe the localized amplitude as that obtained from the unrestricted SD cases, but 
the nodal structure is essentially the same as that obtained from the restricted SD. In order 
to obtain a variationally lower energy we have tried a restricted MRCI (multi reference CI) 
using orbitals obtained from the present CASSCF calculation, but QMC calculations using 
the MRCI trial/guiding wavefunction could not give better results than when using the 
UB3LYP-SD one. 



C. Unrestricted multi determinant calculations 



Several advanced MD implementations such as CASSCF are available at the restricted 
level, but we could not obtain better results than UB3LYP-SD. We therefore tried unre- 
stricted CI (UCI) methods, for which Gaussian03 [58( was used to provide UCISD (UCI 
singles and doubles). Using a UHF reference, the UCI expansion gives 7,521,823 terms, all 
of which can not be taken into account in a QMC calculation. We truncated this expression 
into 35 determinants, removing those terms with coefficients |c$| < 0.01. The expansion 
coefficients were optimized further by VMC, first by weight-limitted variance minimizations 



591 ]. followed by energy minimizations. The results at the experimental R e are shown in 



Table IVHI At the SCF level, UCISD achieved a lower energy than its initial guess UHF, 
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implying that the coefficients are well optimized by CI-SCF. UCISD gives a higher energy 
than UB3LYP at the SCF level, which does not matter, because there is no variational 
relation between them. Using the UCISD trial wavefunction, we achieved a lower VMC 
energy than when using the UB3LYP one, indicating that our VMC optimization of the CI 
coefficients was successful. At DMC, however, UCISD turned out to give a worse result than 
UB3LYP. We also tried another choice of expansion: 67 excited configurations, in which the 
active space was HOMO±6 and the single and double excitations of occupied orbitals were 
restricted to virtual orbitals with the same symmetry. This choice gives, however, a worse 
result than UB3LYP-SD even at the VMC level [-172.744(2) hartree]. 

The above CI treatments could not give any variationally better trial node than UB3LYP- 
SD. The easiest way to go beyond those treatments would be to add UCI expansions to 
UB3LYP-SD because it is the best starting point. By considering the UB3LYP orbital 
symmetry near the HOMO level, we made two different sizes of CI expansions: The first one 
took into account only 3 virtual orbitals above the LUMO, including only a and 7r symmetries 
[for which we refer it as UCISD (UB3LYP+3)], and the second was a larger one with 10 
virtual orbitals in which a, ir, and 5 symmetries were included [UCISD(UB3LYP+10)]. 
In both cases we considered only such excited configurations between the orbitals with 
the same symmetry, resulting in around 50 and 650 determinants for UCISD (UB3LYP+3) 
and UCISD(UB3LYP+10), respectively (the numbers of determinants vary a little amount 
depending on R). 

The results are shown in Fig. |6] (a) and (b). As seen in Fig. |6] (a), VMC using the 
UCISD (UB3LYP+3) trial wavefunction gives a better result than when using the UB3LYP- 
SD one, because the former includes more variational degrees of freedom to be optimized. For 
UCISD(UB3LYP+10), however, we could not get a satisfactory optimization, giving a higher 
energy than the initial UB3LYP-SD calculation. Nevertheless the DMC calculation with the 
UCISD(UB3LYP+10) node gives a slightly lower energy than that with the UB3LYP node 
(Fig. El (b)). It is apparent that the UCISD(UB3LYP+3) node give a lower DMC energy 
than the UB3LYP-SD node. Focussing on UCISD (UB3LYP+3), we further introduced the 
backflow transformation, getting the best binding curve beyond SD, as shown in Fig. [7] 
(a) and (b). Similar to the SD cases, the DMC projection gets R e shorter and D e deeper 
than VMC. Though UCISD (UB3LYP+3) gives a variationally better result than any SD 
treatment, it could not hardly improve the binding nature, i.e., it overestimates R e and 
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underestimates D e , compared with the experimental values. As seen in Fig. 13(a) and (b), 
the backflow does not improves these as well. 

IV. CONCLUDING REMARKS 

We studied the binding curve of the ground state C12 dimer using the fixed node DMC 
method. Various different types of nodal structures were compared based on the variational 
principle with respect to the node of the DMC guiding function. We tested several choices of 
XC functionals with or without spin restriction on the orbital functions composing the many- 
body wavefunction. We also tried computationally expensive choices of UCI expansions with 
backflow transformation. 

Within the SD treatment, UB3LYP turns out to give the variationally best trial node. 
Except HF, the unrestricted nodes are found to be better than the restricted ones. Any 
choice gives binding curves with a energy minimum, but for unrestricted trial nodes they 
end up with a much larger R e and smaller D e in the QMC final results, compared with 
the experimental values. Though some unrestricted DFT-SCF calculations, such as ULDA 
and UGGA, reproduce a proper R e , it is found that the QMC calculations with these trial 
nodes overestimate R e . The restricted nodes recover fairly well R e even at QMC, but they 
give a higher energy than the unrestricted ones. At the experimental R e , we could not get 
a stable molecular energy lower than twice the atomic energy at the QMC level, although 
some DFT-SCF calculations did reproduce a bound state. 

We also examined whether the binding curve could be improved by different pseudopo- 
tentials, comparing the BFD potential with the LN one. Both potentials give almost the 
same binding curve, which justifies our choice of pseudopotentials. The time step bias is 
confirmed to be kept within the error bar considered in this study. The backflow transforma- 
tion turned out to give no specific improvements on describing the binding nature, although 
it did improve the energy variationally. 

Within the framework of the MD approximation, we first tried the GVB model as a 
compact expansion for the many-body wavefunction because of its plausible physical mean- 
ing. Starting with HF orbitals, we optimized the GVB orbitals and coefficients by a SCF 
procedure, but it could not give a better result than the unrestricted SD in the DMC final 
results, probably because the restricted treatment has a limitation in describing the spin 
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polarized nature in Ci2- We also run DMC using the trial nodes derived from a restricted 
MD expansion with orbitals optimized by CASSCF, but it ends up with a worse QMC re- 
sult than the unrestricted DFT-SD nodes. Though the restricted CASSCF itself gives a 
too large R e similar to the unrestricted DFT calculations, the QMC calculations with the 
CASSCF nodes give properly shorten R e near to the experimental value. Then we decided 
to use the unrestricted CI nodes. First, we chose such excited configurations as those with a 
large coefficient weight and re-optimized the coefficients. VMC with this trial wavefunction 
gives a lower energy than that with the UB3LYP one. On the other hand, the corresponding 
DMC calculation gives a higher energy than DMC using the UB3LYP nodes. Next, we man- 
ually constructed a UCI expansion with UB3LYP orbitals and optimized their coefficients. 
Though the UCI-MD node gives a better DMC result than the UB3LYP-SD one, it gives no 
improvements on the binding natures, i.e., it still overestimates R e and underestimate D e . 
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Appendix A: Jastrow functions 



In CASINO 59j, u, x, and / terms in Eq. flj]) are given in power expansion form as 62 1 
u ( r ij) = ( r ij ~ L u)° xQ(L u - r^) 



x \ a Q + 



T; 
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In QWalk 6JJ, instead, the terms are expanded with basis sets {bj (r)} which vanish at some 
cutoff length , given as 



u (r^) = u (rij) + £ c m • &< 
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(A8) 
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(A10) 



where the upper index (pq) in Eqs. IA9l and lAl0l stands for particle pairs such as (el) or (ee). 
The term uq imposes the electron-electron cusp condition for spin pair aa' with a projection 
operator P acr > and cutoff length r a c & ' . The Poly-Pade type basis bm q ^ with (r) is designed 
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to be cutoff quadratically at &{, . The parameters are generated from a given (3q 

using the following recursions: 

ft = exp(1.6) x (l + /3 ) , (All) 
& = [exp(1.6)] xft_i (fc > 1) . (A12) 

The range of summation for the /-term, (k, m, n), denotes the pairs specified according to a 
given order (M u , M x ). For the present work (M u , M x ) = (3, 3) amounts to 12 terms for the 
expansion. 

For CASINO we used expansion orders of A" u =8, N x =8, N^ N =2, and N^=2, with fixed 
cutoff lengths L u = 5 [a.u.], L x = 4 [a.u.], and Ln = 3[a.u.], respectively. In QWalk we 
choose cutoff lengths to be fixed as 7.5 [a.u.]. 
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FIG. 1: Binding curves obtained from unrestricted (a) SCF, (b) VMC, and (c) DMC calculations. 
In (a), SCF energies are shifted so that each minimum is set to be zero; the minimums are, 
-172.758,-173.208,-173.026,-172.946, and -173.026 hartree for SVWN5, PW91PW91, BLYP, 
B1LYP, B3LYP, respectively. In (b) and (c) error bars are within symbol size. VMC and DMC 
include the backflow transformation. 
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FIG. 2: Binding curves obtained from restricted (a) SCF, (b) VMC, and (c) DMC calculations. 
In (a), SCF energies are shifted so that each minimum is set to be zero; the minimums are, 
-172.753,-173.180,-173.004,-172.794, and -172.900 hartree for SVWN5, PW91PW91, BLYP, 
B1LYP, B3LYP, respectively. In (b) and (c) error bars within symbol size. VMC and DMC include 
the backflow transformation. 
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(a) SCF/UB3LYP/BF 
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FIG. 3: Comparison of binding curves using different pseudopotentials with UB3LYP/BF for (a) 
SCF, (b) VMC, and (c) DMC calculations. Energies are shifted so that each minimum to be zero: 
the values of the minimum are (a) —173.026 hartree for LN and —173.849 hartree for BFD, (b) 
-172.920 hartree for LN and -173.746 hartree for BFD, and (c) -173.029 hartree for LN and 
-173.844 hartree for BFD. 
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FIG. 4: Time step dependence of ground state energy in the DMC calculation with the backflow 
transformation. The fixed node is generated from the UB3LYP-DFT calculation. 
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FIG. 5: (a) SCF, (b) VMC, and (c) DMC binding curves using the CASSCF(2,4), CASSCF(2,7), 
and UB3LYP (trial) wavefunctions. In (a) SCF energies are shifted so that each minimum, 
— 171.540, —171.461, and —173.026 hartree to be zero. In (b) and (c) error bars are within symbol 



size. 
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(a) VMC/UB3LYP and its UCISD/without Backflow 
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FIG. 6: Comparison of QMC binding curves using the UB3LYP-SD and UCISD(UB3LYP+iV) 
(N = 3, 10) trial/guiding wavefunctions for (a) VMC and (b) DMC without Backflow transforma- 
tion. 
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FIG. 7: Comparison of QMC binding curves using the UB3LYP-SD and UCISD(UB3LYP+3) 
trial/guiding wavefunctions for (a) VMC and (b) DMC with the backflow transformation. 
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(a) Effect of Backflow at VMC 
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FIG. 8: Comparison between results with and without the backflow transformation for (a) VMC 
and (b) DMC. 



27 



TABLE I: Exchange-Correlation potentials. 'GC stands for gradient correction. 







Exchange 




Correlation 




Non-local 


Local 


GC 


Local 


GC 


Functional 


Vf F [%] 


y| latcr [%] 


8V X [%] 


F VWN5 [%] 


SV C {%} 


SVWN5 





100 





100 





PW91PW91 





100 


100 PW91 


100 


100PW91 


BLYP 





100 


100 B88 


100 


100 LYP 


B1LYP 


25 


75 


75B88 


100 


100 LYP 


B3LYP 


20 


80 


72B88 


100 


81 LYP 
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TABLE II: Ground state energies [hartree] at the experimental bond length. Underlined values 
highlight those lower than the zero-binding energy in Table \TU\ 





SCF 


VMC 


DMC 


UHF 


-171.701 


-172 


665(6) 


-172.903(2) 


RHF 


-171.097 


-172 


584(4) 


-172.911(2) 


USVWN5 


-172.757 


-172 


765(2) 


-172.976(2) 


RSVWN5 


-172.741 


-172 


704(2) 


-172.963(3) 


UPW91PW91 


-173.201 


-172 


796(2) 


-172.983(2) 


RPW91PW91 


-173.173 


-172 


681(2) 


-172.954(2) 


UBLYP 


-173.022 


-172 


792(2) 


-172.985(2) 


RBLYP 


-173.000 


-172 


719(2) 


172.965(1) 


UB1LYP 


-172.887 


-172 


800(2) 


-172.974(3) 


RB1LYP 


-172.777 


-172 


692(2) 


-172.952(2) 


UB3LYP (without BF) 


-172.974 


-172 


756(2) 


-172.954(2) 


UB3LYP(with BF) 


-172.974 


-172 


819(2) 


-172.985(2) 


RB3LYP 


-172.885 


-172 


710(2) 


-172.963(2) 


ZeroBinding(Best) 




-172 


931(2) 


-173.012(3) 



TABLE III: Zero-binding energies [hartree] for different trial wavefunctions. QMC values are 
evaluated with the backflow transformation. 





SCF 


VMC 


DMC 


SVWN5 


-172.668 


-172.918(2) 


-173.007(1) 


PW91PW91 


-173.163 


-172.919(2) 


-173.007(1) 


BLYP 


-172.972 


-172.925(2) 


-173.008(2) 


B1LYP 


-172.913 


-172.931(2) 


-173.012(3) 


B3LYP 


-172.972 


-172.910(2) 


-173.000(2) 


HF 


-171.875 


-172.921(2) 


-173.002(2) 
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TABLE IV: Comparison of zero-binding energies [hartree] for different pseudopotentials. SCF 
molecular energies are also listed for comparison, designated as SCF (molecule). 



pseudopotentials 


SCF 


SCF (molecule) 


VMC 


DMC 


Lee-Needs (LN) 


-172.972 


-172.974 


-172.910(2) 


-173.000(2) 


Burkatzki-Filippi-Dolg (BFD) 


-173.809 


-173.780 


-173.740(2) 


-173.820(2) 



TABLE V: Symmetries of the UHF natural orbitals near the HOMO (14) and LUMO (15) levels. 



Level 


Symmetry 


20 


<74s 


19 


xz 


18 


yz 


17 


a z 2 


16 


xy 


15 


x 2 — y 2 


14 


x 2 - y 2 


13 


xy 


12 


a z 2 


11 


yz 


10 


xz 


9 


Pis 



30 



TABLE VI: GVB energies [hartree] at experimental bond length. 'GVB(6)opt' stands for that 
with coefficients optimized further by VMC (see text). 



Methods 


SCF 


VMC 


DMC 


UHF 


-171.701 


-172.665(6) 


-172.903(2) 


RHF 


-171.097 


-172.584(4) 


-172.911(2) 


GVB(6) 


-171.624 


-172.698(6) 


-172.861(2) 


GVB(6)opt 




-172.723(2) 


-172.933(2) 



TABLE VII: Comparison between UCISD-MD and UB3LYP-SD at experimental bond length. 
QMC results are those without the backflow transformation. Underlined value highlights that 
reproducing the binding. 



Methods 


SCF 


VMC 


DMC 


UHF 


-171.701 


-172.665(6) 


-172.903(2) 


UB3LYP(w/o BF) 


-172.974 


-172.756(2) 


-172.954(2) 


UCISD(w/o BF) 


-172.599 


-172.772(2) 


-172.926(3) 


ZeroBinding(Best) 




-172.931(2) 


-173.012(3) 
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